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Percolation on two-dimensional small- world networks has been proposed as a model for the spread 
of plant diseases. In this paper we give an analytic solution of this model using a combination of 
generating function methods and high-order series expansion. Our solution gives accurate predic- 
tions for quantities such as the position of the percolation threshold and the typical size of disease 
outbreaks as a function of the density of "shortcuts" in the small-world network. Our results agree 
with scaling hypotheses and numerical simulations for the same model. 



I. INTRODUCTION 

The small-world model has been introduced by Watts 
and Strogatz ||l| as a simple model of a social network — 
a network of friendships or acquaintances between in- 
dividuals, for instance, or a network of physical con- 
tacts between people through which a disease spreads. 
The model consists of a regular lattice, typically a one- 
dimensional lattice with periodic boundary conditions 
although lattices of two or more dimensions have been 
studied as well, with a small number of "shortcut" bonds 
added between randomly chosen pairs of sites, with den- 
sity (j) per bond on the original regular lattice. The 
small- world model captures two specific features observed 
in real- world networks, namely (1) logarithmically short 
distances through the network between most pairs of in- 
dividuals and (2) high network clustering, meaning that 
two individuals are much more likely to be friends with 
one another if they have one or more other friends in 
common. The model turns out to be amenable to treat- 
ment using a variety of techniques drawn from statistical 
physics and has as a result received wide attention in the 
physics community ||, |, |, |, §, 0]. 

The small- world model has some problems however. In 
particular, it is built on a low-dimensional regular lattice, 
and there is little justification to be found in empirical 
studies of social networks for such an underlying struc- 
ture. As recently pointed out by Warren et al. how- 
ever, there is one case in which the small-world model 
may be a fairly accurate representation of a real-world 
situation, and that is in the spread of plant diseases. 
Plant diseases spread through physical contacts between 
plants — immediate contagion, insect vectors, wind, and 
so forth — and these contacts form a "social network" 
among the plants in question. However, plants are ses- 
sile, and confined by and large to the two-dimensional 
plane of the Earth's surface. Disease spread as a re- 
sult of short-range contact between plants is thus proba- 



bly well represented as a transmission process on a sim- 
ple two-dimensional lattice, and disease spread in which 
some portion of transmission is due to longer-range vec- 
tors such as wind or insects may be well represented 
by a small-world model built upon an underlying two- 
dimensional lattice. Because of this, as well as because 
of inherent mathematical interest, a number of recent 
papers have focused on two-dimensional small-world net- 
works |, 10. 

In this paper we study bond percolation on 
two-dimensional small-world networks. Bond per- 
colation is equivalent to the standard suscepti- 
ble/infectious/recovered (SIR) model of disease 
spread 11 , in which all individuals are initially 
susceptible to the disease, become infected (and hence 
infectious) with some probability per unit time if one of 
their neighbors is infectious, and recover again, becoming 
uninfectious and also immune, after a certain time in 
the infectious state. The equivalence of the SIR model 
to percolation is straightforward, with the percolation 
threshold mapping to the epidemic threshold of the 
disease in the SIR model, and cluster sizes mapping 
to the sizes of disease outbreaks which start with a 
single disease carrier. Using a combination of an exact 
generating function method with a high-order series 
expansion, we derive approximate analytic results for 
the position of the threshold and the mean outbreak and 
epidemic sizes as a function of the density of shortcuts 
in the 2D small-world network and the percolation 
probability, which is equivalent to disease transmission 
probability. As we demonstrate, our results arc in 
excellent agreement with those from other studies using 
different methods, as well as with our own numerical 
simulations. 
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II. GENERATING FUNCTION FORMALISM 



self-consistently as [O 



We study the two-dimensional small- world model built 
on the square lattice. The model is depicted in Fig. ||. 
We develop our generating function formalism for the 
general case of a d-dimensional square/cubic/hypercubic 
underlying lattice first, narrowing our scope to the two- 
dimensional case in Section III where we describe our 



series expansion calculations. For a d-dimensional lattice 
with bonds along the principal axes out to distance k , 
the underlying lattice has dL'^k bonds on it, where L 
is the system dimension, and shortcuts are added with 
probability cj) per underlying bond, for a total of dL'^kcj) 
shortcuts. Then all bonds, including the shortcut bonds, 
are occupied with probability p, or not with probability 
1 — p, and we construct the percolation clusters of sites 
connected by the occupied bonds. 

The generating function part of our calculation follows 
the method of Moore and Newman jl2j. We define a 
probability generating function H{z) thus: 



iI(z) = ^P(n)^", 



(1) 



where P{n) is the probability that a randomly chosen 
site in our small-world network belongs to a connected 
cluster of n sites other than the system-spanning cluster. 
Note that if the probability distribution P{n) is prop- 
erly normalized, then H{1) = 1 below the transition and 
H{1) = 1 — S above the transition where S is the fraction 
of the system occupied by the system-spanning cluster. 

We also define Po{n) to be the probability that a ran- 
domly chosen site belongs to a cluster of n sites on the 
underlying lattice. The complete cluster on the small- 
world network is composed of a set of such underlying 
clusters, joined together by occupied shortcut bonds. If 
we denote by P(m\n) the probability that an underly- 
ing cluster of n sites has exactly m shortcuts emanating 
from it, then the generating function H(z) can be written 



(a) 





FIG. 1: (a) A two-dimensional small-world network built 
upon a square lattice with connection range k = 1. (b) When 
fc > 1 the underlying lattice contains bonds beyond nearest- 
neighbor bonds along each principal axis out to range k, as 
shown here for A; = 3. 



Hiz) 



^PoHz"5]P(m|n)[i/(z)]^ 



71=1 



(2) 



(The derivation of this equation assumes that all clusters 
other than the percolating cluster, if there is one, contain 
no closed loops other than those on the underlying lattice, 
i.e., there are no loops involving shortcut bonds. This is 
only strictly true in the limit of infinite system size, and 
hence our results will only be exact in this limit.) 

There are a total of 2dL'^k(j>p ends of occupied short- 
cuts in the model, and since all ends are uniformly dis- 
tributed over the lattice the probability that any end 
lands in a given cluster of size n is just n/L'^ for large L. 
Thus P{m\n) is given by the binomial distribution 



P{m\n) 



l2dL'^k(f)p\^ 


n 


m 






1? 





(3) 



Substituting into Eq. (^) and performing the sum over 
TO, this gives 



71 



l + {H{z)-l)-^ 



-1 2dL''ktpp 



(4) 



where the last equality holds in the limit of large L. If 
we define the additional generating function 



Ho{z)^J2Poin)^ 



(5) 



which is the probability generating function for the sizes 
of clusters for ordinary bond percolation on the underly- 
ing lattice, then Eq. (||) can be written in the form 



H{z)^Ho{ze^'""^P^"^'^-'^). 



(6) 



This gives a self-consistency condition from which we can 
evaluate H{z), and hence we can evaluate the probabili- 
ties P{n) for cluster sizes in the small-world model. 

In fact, it is rarely possible to solve Eq. (||) for H{z) in 
closed form (although evaluation by numerical iteration 
is often feasible), but we can derive closed- form expres- 
sions for other quantities of interest. In particular, the 
average size of the cluster to which a randomly chosen 
site belongs is given by 



in) 



iP{n) = H'{1) = i?o(l)[l + 2dk(PpH'{l)], 



or, rearrangmg. 



1 - 2dk(l)pH'^{iy 



(7) 



(8) 
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This quantity diverges when 



2dA;#i?^(l) 



or equivalently when 



2dkpH'^{iy 



(9) 



(10) 



and this point marks the phase transition at which a gi- 
ant cluster scahng as the size of the entire system first 
forms. Another way of looking at this result is to note 
that Hq{1) = (no), the average cluster size on the under- 
lying lattice. Thus percolation takes place when 



2dk(j)p = 



1 



(11) 



The quantity on the left-hand side of this equation is the 
average density of the ends of occupied shortcut bonds 
on the lattice (see Section IV), and thus percolation takes 
place when there is exactly one end of an occupied short- 
cut bond per cluster on the underlying lattice. This is 
reminiscent of the phase transition in an Erdos-Renyi 
random graph, which occurs at the point where each ver- 
tex in the graph is attached to exactly one edge . 

Above the phase transition S = 1 — H{1) is the size of 
the giant cluster. Setting z = 1 in Eq. (||) we find that S 
is a solution of 



S=l-H, 



o(c 



'2dk4>pS\ 



(12) 



which can be solved by numerical iteration starting from 
a suitable initial value of S. An expression similar to (|^) 
for the average size of non-percolating clusters above the 
transition can also be derived. See, for example, Ref. [ij. 



III. SERIES EXPANSIONS 



In the one-dimensional case studied in Ref. 12, the cal- 



culation of the generating function Hq{z) is trivial — it 
is equivalent to solving the problem of bond percolation 
in one dimension. In the present case however we are 
interested primarily in the two-dimensional small-world 
model, and calculating iJp is much harder; no exact so- 
lution has ever been given for the distribution of cluster 
sizes for bond percolation on the square lattice. Instead, 
therefore, we turn to series expansion to calculate i/o 
approximately. 

Ho{z) for the two-dimensional case can be written as 



Hoiz) 



stn 



(13) 



where gstn is the number of different possible clusters 
on a square lattice which have s occupied bonds, t un- 
occupied bonds around their perimeter, and n sites. If 
we can calculate gstn up to some finite order then we 
can calculate an approximation to Hq{z) also. Unfortu- 
nately, because gstn depends separately on three different 



indices, it is prohibitively memory-intensive to calculate 
on a computer up to high order. We note however that 
we only need Hq as a function of p and z, and not of 1 — p 
separately, so we can collect terms in p and rewrite the 
generating function as 



H,{z)^Y.P^ 



(14) 



where the quantities Qm{z) are finite polynomials in z 
which are, it's not hard to show, of order z™+^. Calcu- 
lating these polynomials is considerably more economical 
than calculating the entire set of gstn- We have calcu- 
lated them up to order to = 31 using the finite-lattice 
method , in which a generating function for the infi- 
nite lattice is built up by combining generating functions 
for the same problem on finite lattices. The finite lat- 
tices used in this case were rectangles of h x I sites and 
the quantity we consider is the fundamental generating 
function for the cluster density 



G(z,p) = ^zV(l-p)* 



5str: 



(15) 



which, it can be shown, is given by the linear combination 
G{z,p) = ^WhiGhi{z,p), (16) 

hi 

where Whi are constant weights that are independent of 
both z and p, and Ghi{z,p) is the generating function for 
connected clusters (bond animals) which span an h x I 
rectangle both from left to right and from top to bottom. 

Due to the symmetry of the square lattice, the weight 
factor Whi is simply 



Whl = 



for I < h 
for I = h 
for ; > h. 



(17) 



The individual generating functions Ghi{z,p) for the 
finite lattices are calculated using a transfer matrix 
method with generating functions for all rectangles of 
a given height h being evaluated in a single calculation 
The algorithm we use is based on that of Conway |16 
with enhancements similar to those used by Jensen 1 17 
for the enumeration of site animals on the square lat- 
tice |l^. Since clusters spanning a rectangle of h x I 
sites contain at least h + I ~ 2 bonds, we must calculate 
Ghi{z,p) for all h< {m+l)/2 and /i<Z<TO-/i-|-2in 
order to derive a series expansion for G{z,p) correct to 
order m in p. For the order m = 31 calculation described 
here the maximal value of h required was 16. 

Once G{z,p) is calculated, the polynomials Qm{z) are 
easily extracted by collecting terms in p. (Alternatively, 
one could write Ho{z) = zdG/dz, although doing so of- 
fers no operational advantage in the present case.) In Ta- 
ble I we list the values of Qmiz) for to up to 10; the com- 
plete set of polynomials up to to = 31 is available from 
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m 


Qmiz) 





Z 


1 


4z^ ~ 4z 


2 


18z'^ — 242^ + 6z 


3 


882* — 1442'' + 602^ — 42 


4 


4352^ — 8602^ + 5042^ - 802^ + z 


5 


21842*' — 50202^ + 37842* - 10082^ + 602^ 


R 
\J 


1 1 ni Sz^ — 98Q'^9y^ -1- 9fi'i'=iny^ — qS79r'* -1- 1 9fin7^ — 947-^ 

_L_LL'X02 jIOUO^Z ^ ^\j00\JZ — yO ( ±^UU2 — Zj^/C 


7 


558882* - 1646682'^ + 1779722*' - 851002^ + XmVlz^ - 10082=* + 42^ 


8 


2842292^ - 9288402** + 11536982^ - 6738362^ + 1841252^ - 198802^' + 5042^ 


9 


14488002"' - 51971762^ + 72914882® - 50303122'^ + 17544242^ - 2833202^ + 162402'' - 1442^ 


10 


73962902" - 288901602"' + 451559522^ - 359267202* + 152788722'^ - 33230882^ + 3179402^ - 91042* + 182^ 



TABLE I: The values of the polynomials Qm(z') up to m = 10. 



the authors on request. We notice that since -ffo(l) — 1 
for all p < Pci as it must given that the probability dis- 
tribution it generates is properly normalized, it must be 
the case that 



.(1) 



1 for TO = 
for TO > 1. 



(18) 



It can easily be verified that this is true for the orders 
given in Table |. This implies that ffo will be correctly 
normalized even if we truncate its series at finite order 
in p, as we do here. This makes our calculations a little 
easier. 

In order to calculate the average cluster size and 
position of the phase transition ( [lO|) in our small-world 
model, we need to evaluate the quantity Hq{1), which is 
given by 



m=0 



(19) 



The quantities Q'„-^{1) are just numbers — their values up 
to TO = 31 are given in Table |^ — so that this expression 
is a simple power series in p. If we make use of our results 
for TO up to 31 to evaluate this quantity directly, we can 
calculate the behavior of the 2D small- world model using 
the results of Section ||. However we can do better than 
this. 

Since Hq{1) is the average size (no) of a cluster in or- 
dinary bond percolation on the square lattice, we know 
that it must diverge at pc = ^, and that it does so as 
{pc — p)~'^ , where 7 is the mean cluster-size exponent for 
two-dimensional percolation which is equal to y|. With 
this information we can construct a Pade approximant 
toi/^(l) [|0|,|l|. Writing 



Pc-P 



Pc 



(20) 



where A{p) is assumed analytic near pc, we construct a 
Pade approximant to the series for 



A{p)^ 



Pc-P 
Pc , 



(21) 



using our series for Hq(1). Then we use this approximant 
in Eq. ( pO| ) to give an expression for Hq{1) which agrees 
with our series expansion result to all available orders, 
and has a divergence of the expected kind at p = i. 
As is typically the case with Pade approximants, the 
best approximations are achieved with the highest order 
symmetric or near-symmetric approximants and using all 
available orders in our series expansion, we find the best 
results using a [15, 15] approximant to A{p) in Eq. (pi]). 
In Fig. ||we show the resulting estimate for (no) = Hq{1) 
(dotted line) as a function of p against numerical results 
for the average cluster size on an ordinary square lattice 
(squares). As the figure shows, the agreement is excel- 
lent. 

Substituting our Pade approximant expression for 
Hq{1) into Eq. (||) we can now calculate average clus- 
ter size for the small- world model for any value of cf), and 



m 


Q™(1) 


m 


Q™(1) 





1 


1 


4 


2 


12 


3 


36 


4 


88 


5 


236 


6 


528 


7 


1392 


8 


2828 


9 


7608 


10 


14312 


11 


39348 


12 


69704 


13 


197620 


14 


318232 


15 


1013424 


16 


1278912 


17 


5362680 


18 


4418884 


19 


28221636 


20 


11543548 


21 


152533600 


22 


-20880672 


23 


903135760 


24 


-705437704 


25 


5680639336 


26 


-7577181144 


27 


37205966052 


28 


-66485042424 


29 


253460708032 


30 


-534464876516 


31 


1767651092388 



TABLE II: The derivatives g^(l) for all orders up to m = 
31, which are also the coefficients of the series for the mean 
cluster size — see Eq. ([ig|). Note that although the values for 
Q'mW appear initially to be positive and increasing roughly 
exponentially, this rule does not hold in general. The first 
negative value is at m = 22, and the signs of Q'mi^) appear 
to alternate for m > 22. 



5 




density of shortcuts (]) 

occupation probability p 



FIG. 2: The average of size in sites of the cluster to which a 
randomly chosen site belongs for bond percolation on a two- 
dimensional small- world network with k = 1. The circles are 
simulation results for systems of 1024 x 1024 sites, calculated 
using the fast algorithm of Newman and Ziff |19|, and the 
solid lines are the analytic result, Eq. (^. From left to right, 
the values of for each of the lines are 1.0, 0.5, 0.2, 0.1, 
0.05, and 0.02. As (j> diminishes the lines asymptote to the 
normal square lattice form which is indicated by the square 
symbols (simulation results) and the dotted line (series ex- 
pansion/Pade approximant result). 



from Eq. (^0|) we can calculate the position of the phase 
transition. In Fig. ^ we show the results for (n) as a 
function of p along with numerical results for the same 
quantity from simulations of the model. In Fig. ^ we 
show the results for pc plotted against numerical calcula- 
tions In both cases, the agreement between analytic 
and numerical results is excellent. In the inset of Fig. ^ 
we show the results for pc on logarithmic scales, along 
with the value calculated by using the series expansion 
for H'q{1) directly in Eq. (|l^ (dotted hue). As the figure 
shows, the Fade approximate continues to be accurate to 
very low values of 0, where the direct series expansion 
fails. 



IV. SCALING FORMS 

As Ozana has pointed out, there are two com- 
peting length-scales present in percolation models on 
small- world networks. One is the characteristic length 
^ of the small-world model itself which is given by ^ = 
\/{2(j)kdY^'^^ where d is the dimension of the underly- 
ing lattice (2 in the present case). This length is the 
typical linear dimension of the volume on the underlying 
lattice which contains the end of one shortcut on aver- 
age. In other words = 2(j>kd is the density of the 
ends of shortcuts on the lattice. (In fact, one normally 
leaves the factor of 2 out of the definition of the char- 



FIG. 3: The position of the percolation transition for the 2D 
small-world network as a function of the density (p of short- 
cuts. The points are simulation results, again using the al- 
gorithm of Ref. |l^, and the lines are our analytic calculation 
using a Fade approximant. For the simulations, the value of 
Pc was taken to be the point at which the size of the largest 
cluster in the system has maximal gradient as a function of p. 
The slight difference between the numerical and analytic re- 
sults appears to be a systematic error in the estimation of Pc 
from the numerical results (see Ref. ^ for a discussion of this 
point). Inset: the same comparison on logarithmic scales. In 
this case, the vertical axis measures | —Pc, which should go to 
zero with (j> according to i — pc ~ 0^''^ where 7 = ||- (This 
can be deduced from Eq. (p^, and is also shown in Ref. ^.) 
Again the solid line is the Fade approximant calculation, while 
the dotted line represents the value of pc calculated directly 
from the series expansion without using a Fade approximant. 

acteristic length, but we include it since it makes the 
resulting formulas somewhat neater in our case.) In the 
current percolation model, only a fraction p of the short- 
cuts are occupied and thereby contribute to the behavior 
of the model — the unoccupied shortcuts can be ignored. 
Thus the appropriate characteristic length in our case is 
derived from the density of ends of occupied shortcuts, 
which was discussed in Section ||. The correct expression 
is 



(2p#rf)i/'i- ^ ^ 

The other length scale is the correlation length or typ- 
ical cluster dimension for the percolation clusters on the 
underlying lattice. Normally the latter is also denoted 
but to avoid confusion we follow the notation of Ref. 
and here denote it Ozana has given a finite-size scaling 
theory for percolation on small- world networks which ad- 
dresses the interaction of each of these length-scales with 
the lattice dimension L. Our analytic calculations how- 
ever treat the case of L — > oo, for which a simpler scaling 
theory applies. In this case, the only dimensionless com- 
bination of lengths is the ratio C/Cy any observable 
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quantity Q must satisfy a scaling relation of the form 

Q ^ c'^r'^'/ic/c)- (23) 

where a and [3 are scaling exponents and ,f{x) is a uni- 
versal scaling function. This form applies when we are 
in the region where both ^ and C are much greater than 
the lattice constant, i.e., when shortcut density is low 
(the scaling region of the small-world model) and when 
we are close to the percolation transition on the normal 
square lattice. 

We can rewrite Eq. (|2^) in a simpler form by making 
use of Eq. (|2^) and the fact that the typical cluster di- 
mension C, on the underlying lattice is related to typical 
cluster volume by = {no) (assuming compact clusters). 
This then implies that 



transition in powers of 5, to give 



Q - {noT{p4>kdfF{2p^kd{no)), 



(24) 



where F{x) is another universal scaling function. 

Consider for example the average cluster size (n) for 
the small- world model. Since we know this becomes equal 
to its normal square-lattice value (no) when </() = 0, we 
can immediately assume a — 1, (3 — and 



F{2p4)kd{no)). 



(25) 



Thus, a plot of (n) / (no) against the scaling variable 
X = 2p(j)kd{nQ) should yield a data collapse whose form 
follows the scaling function F{x). In fact there is no need 
to make a scaling plot in this simple case. Comparison of 
Eq. ( |25| ) with Eq. (||), bearing in mind that {uq) = Hq{1), 
reveals that (n) does indeed follow the expected scaling 
form with 



F{x) = 



1 



1 



(26) 



The point x = 1, which is also the point at which the 
two length-scales are equal ^ = Ci thus represents the 
percolation transition in this case. (This observation is 
equivalent to Eq. (pT|).) 

A slightly less trivial example of the scaling form (|2^) is 
the scaling of the size S of the giant percolation cluster, 
Eq. (p^. To deduce the leading terms in the scaling 
relation for we expand (f2|) close to the percolation 



S = xS 



[jm+Hiii)] 



x'S' + OiS-"). 



(27) 



Rearranging and keeping terms to leading order, we find 



S: 



[i/^(l)-fi/^'(l)] 



(28) 



If there is only one correlation length for percolation on 
the underlying lattice then {na)"^ / {n^) is homogeneous in 
it and hence constant in the critical region. Thus S scales 
as 5 ~ {x — l)/x'^ close to the transition, with the leading 
constant being zero below the transition and 0(1) above 
it. 

V. CONCLUSIONS 

We have presented analytic results for bond percola- 
tion in the two-dimensional small-world network model, 
which has been proposed as a simple model of the spread 
of plant diseases. Using a combination of generating func- 
tion methods and series expansion, we have derived ap- 
proximate but highly accurate expressions for quantities 
such as the position of the percolation transition in the 
model, the typical size of non-percolating clusters, and 
the typical size of the percolating cluster. Our results 
are in excellent agreement with numerical simulations of 
the model. By judicious use of Pade approximants for 
the series expansions, the results can even be extended 
to very low shortcut densities, where a simple series ex- 
pansion fails. The results are also in good agreement 
with the expected scaling forms for the model. 
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